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I.  Introduction 

This  paper  is  concerned  with  a  multichannel  signal  detection  problem  frequently  encountered 
in  phased-array  radars  and  many  other  applications.  With  extra  spatial  information  provided  by 
multiple  sensors,  higher  performance  of  signal  detection  can  be  achieved  (than  a  single-sensor 
system),  especially  in  detection  of  signals  buried  in  a  background  of  directional  jammers  and 
space-time  correlated  clutter.  A  widely  explored  technology  for  multichannel  signal  detection 
is  space-time  adaptive  processing  (STAP)  [1],  first  proposed  by  Brennan,  Reed  and  Mallett 
[2],  Most  STAP-based  methods,  such  as  the  adaptive  matched  filter  (AMF)  [3]  and  Kelly’s 
generalized  likelihood  ratio  test  (GLRT)  [4],  need  to  invert  a  large  space-time  covariance  matrix. 
These  methods  require  not  only  a  large  number  of  independent,  identically  distributed,  signal- 
free  training  data  to  estimate  the  matrix,  but  they  also  incur  a  high  computational  cost  for  matrix 
estimation  and  inversion. 

A  parametric  STAP  detector  based  on  a  multichannel  autoregressive  (AR)  disturbance  model 
has  been  proposed  for  airborne  radar  applications  [5],  [6]  to  reduce  both  the  training  data 
requirement  and  computation  load.  This  method  is  called  the  parametric  adaptive  matched  filter 
(PAMF)  [6].  While  the  PAMF  detector  has  been  found  to  yield  exceptional  performance  with 
significantly  reduced  training  and  computational  requirements  when  compared  with  fully  adaptive 
STAP  detectors,  the  connections  between  the  PAMF  and  other  reduced-dimensional  or  partially 
adaptive  STAP  detectors  [1],  which  have  similar  benefits  in  training  and  complexity,  remains 
unclear. 

This  paper  aims  to  provide  some  insights  into  this  problem  by  employing  the  conjugate- 
gradient  (CG)  method  to  solve  the  linear  prediction  problem  underlying  the  temporal  whitening 
phase  of  the  PAMF  detector.  Our  choice  of  the  CG  method  is  motivated  by  several  factors.  First, 
as  will  be  shown,  the  CG  algorithm  naturally  leads  to  a  subspace  interpretation  of  the  PAMF 
detector,  and  offers  a  connection  to  the  other  reduced-rank  STAP  detectors.  Second,  the  CG 
method  is  a  computational  efficient  algorithm  to  solve  the  linear  prediction  problem  underlying 
the  PAMF  detector.  In  particular,  for  airborne  radar  applications,  due  to  an  inherent  structure  of 
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the  disturbance  covariance  matrix,  the  CG  algorithm  can  usually  achieve  convergence  using  only 
a  few  iterations,  thus  providing  significant  computational  saving.  Third,  since  the  disturbance 
covariance  matrix  has  a  block-Toeplitz  matrix  structure,  preconditioning  methods(e.g.,  [7],  [8], 
[9])  can  be  employed,  which  are  very  effective  in  further  speeding  up  the  convergence  rate  in 
CG  iterations,  while  adding  up  only  a  modest  computational  overhead  per  iteration  (due  to  the 
block-Toeplitz  structure).  Finally,  as  a  by-product,  we  show  that  the  CG  algorithm  also  yields  a 
new  and  computationally  efficient  AR  model  order  selection  method  that  can  be  integrated  with 
the  CG  iterations. 

The  remainder  of  this  paper  is  organized  as  follows.  The  signal  detection  problem  is  introduced 
in  Section  II.  A  brief  review  of  the  matched  filter  (MF)  and  parametric  matched  filter  (PMF) 
detectors  is  provided  in  Section  III.  In  Section  IV,  the  CG  versions  of  MF  (CG-MF)  and 
PMF  (CG-PMF)  and  a  CG-based  model  order  selection  method  are  proposed.  The  convergence 
rate  of  CG  in  airborne  radar  applications,  along  with  a  preconditioned  CG-PMF  (PCG-PMF) 
detector,  are  also  discussed  there.  In  Section  V,  we  consider  the  adaptive  case  and  present  a  new 
model  order-selection  CG-PAMF  (OSCG-PAMF)  detector,  when  both  the  AR  model  order  and 
coefficients  are  unknown.  The  performance  of  the  proposed  class  of  CG-PMF  and  CG-PAMF 
detectors  is  illustrated  by  numerical  results  in  Section  VI.  Finally  conclusions  are  summarized 
in  Section  VII. 

Vectors  and  matrices  are  denoted  by  boldface  lower-case  and  upper-case  letters,  respectively. 
Transpose,  complex  conjugate  and  complex  conjugate  transpose  are  respectively  represented  by 
(•)T,  (•)*  and  (■)rr ■  C  and  R  denote  the  complex  and  real  number  fields.  CAf((J,.  R)  denotes  an 
additive  multivariate  Gaussian  random  variable  with  mean  vector  fi  and  covariance  matrix  R. 

II.  Data  Model 

Consider  a  received  J-channel  sequence  {x(n)|n  =  1, 2,  •  •  •  ,  N}  corrupted  by  a  space-time 
correlated  disturbance  random  process  c(n).  The  detection  problem  involves  the  following  binary 
hypotheses: 

H0  :  x(n)  =  c(n) 

Hi  :  x(n)  =  as(n)  +  c (n)  (1) 
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where  s (n)  is  a  known  J-channel  signal  and  a  is  its  deterministic  but  unknown  complex 
amplitude.  All  vectors  in  (1)  are  J  x  1  vectors.  For  convenience  of  later  discussions,  de¬ 
fine  the  following  vectors  in  descending  order:  s  =  [sT(N),sT(N  —  l),---  ,  sT(l)]r,  c  = 
[cT (N) ,  cT (N  —  1),  •  ■  ■  ,  cT(l)]T,  x  =  [xT(iV),  x7  (iV  —  1),  ■  ■  ■  ,  xr(l)]T.  It  is  standard  to  assume 
that  the  disturbance  c  is  a  Gaussian  random  vector  with  zero-mean  and  space-time  covariance 
matrix  Rc  E  CJNxJN,  while  the  signal  vector  s (n)  is  deterministic  (Swerling  0  target).  Based 
on  these  assumptions,  x  ~  CJ\T(as,  Rc),  where  a  —  0  under  H0  and  a  /  0  under  //, . 

In  STAP,  the  signal  s  is  known  as  the  space-time  steering  vector: 

s  =  St  ®  ss  (2) 

where  s*  and  ss  denote  the  temporal  steering  vector  and  spatial  steering  vector,  respectively, 
and  ®  denotes  the  Kronecker  product.  For  a  side-looking  uniform  linear  array  (ULA),  we  have 
st  =  (1  /\fN)\el2n(N~l^d  ,■  ■  ■  .  e';27rA,  1] 7  with  a  normalized  Doppler  frequency  f,i  and  ss  = 
(1  /\/j)[el2’r*J“1^s,  •  ■  •  ,  ,  1]T  with  a  normalized  spatial  frequency  fs.  Practically,  the  true 

disturbance  covariance  matrix  Rc  is  unknown,  and  often  an  estimate  can  be  obtained  from  the 
secondary  data: 

1  K 

Rc  =  CkC%  (3) 

k= 1 

where  c*,,  k  —  1, 2  •  ■  •  ,  K,  denote  the  secondary  data  vectors  assumed  to  be  signal-free.  Accord¬ 
ing  to  the  well-known  “RMB”  rule  [10],  we  need  K  >  2JN  —  3  so  that  the  average  output 
signal-to-interference-plus-noise  ratio  (SINR)  loss  caused  by  covariance  estimation  error  is  less 
than  3  dB.  Detectors  with  an  estimated  covariance  matrix  are  often  called  adaptive  methods. 


III.  MF  AND  PMF 


Assuming  a  known  Rc,  the  matched  filter  (MF)  is  obtained  by  maximizing  the  output  SINR 
of  a  linear  receiver  or  the  generalized  likelihood  ratio  (GLR).  The  test  is  given  by  (e.g.,  [3]): 


s"Rr  x 


s^R:^  H0 


<  Vmf 


(4) 


where  //MF  is  the  threshold  of  the  MF.  Equation  (4)  is  the  well  known  matched  subspace  detector 
for  a  rank-1  signal  in  colored  noise.  Consequently,  it  offers  unbeatable  performance  for  the 
detection  problem  considered  in  equation  (1). 
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For  ease  of  exposition,  the  MF  can  also  be  represented  by  using  a  structure  of  temporal 
whitening  cascaded  with  spatial  whitening  arising  from  a  block  LDU  decomposition  of  the 
disturbance  covariance  matrix  [6] .  This  form  of  MF  is  given  by 

|(Q-1/2L-1s)ff(Q-1/2£)|2  _  |s^|2Hi 

(Q-1/2L-1S)^(Q-1/2L-1S)  “  sHs  l„r/MF  (5) 

where  Q  G  CJNxJN  is  a  block-diagonal  matrix  with  Flermitian  matrices  Q(n),  n  —  1,  2,  •  •  •  ,N, 
along  the  main  block  diagonal,  and  L  G  C'/A  y'/  V  is  a  lower  block-triangular  matrix  with  J  x 
J  identity  matrices  along  the  main  block  diagonal.  Both  L  and  Q  come  from  a  block  LDU 
decomposition  of  the  disturbance  covariance  matrix  Rc  =  LQL;/.  Finally, 

(n-l) 

e(n)  =  x(n)  -  J]  Af  (p)x(n  -  p)  (6) 

p= i 


s(n) 


is(n)  —  Q  1/,2(n)e(n) 


Q~1/2(") 


(n-l) 

s(n)  -  J]  Af?(p)8(n-p) 

p= i 


(7) 

(8) 


where  A^{p)  G  CJxJ  is  a  block  element  of  L_1  located  at  the  (n  —  p)-th  block  column  and  the 
n-th  block  row.  Due  to  the  fact  that  there  is  no  performance  penalty  for  the  pre whitening  of  the 
interference  [11,  Ch.6],  (5)  is  equivalent  to  (4). 

If  the  disturbance  c(n)  is  stationary  in  time,  the  MF  can  be  simplified.  A  parametric  matched 
filter  (PMF)  was  introduced  in  [6]  by  modeling  the  disturbance  as  a  stationary  P-th-order 
multichannel  autoregressive  (AR)  process.  Specifically, 

p 

c(n)  =  AH(p)c(n  —  p)  +  £p(n)  (9) 

p=  i 


where  A H(p),p  —  1, 2,  ■  •  ■  ,  P,  is  the  p- th  AR  matrix  coefficient  of  linear  prediction,  and  £p(n ) 


is  the  temporally  white  noise  with  a  spatial  covariance  matrix  Qp.  The  PMF  test  is  given  by  [6] 

■\1V 

vn=P+ 1  ’ 


En=p+1Sp(n)i/p(w)|2gi 

^n=P+lSp(n)Mn) 


%  Vpmf 


where  r^p(n)  =  Q p1^2£p(n)  and 

sp(n)  =  Qp1/2 


s (n)  —  AH(p)s(n  —  p) 

p= i 


(10) 


(11) 


for  n  =  P  +  1,  -  •  ■  ,  A.  In  practice,  the  model  order  P  and  the  AR  coefficients  {A(p)\  are 
unknown  and  hence  estimated  from  the  secondary  data  and/or  primary  data.  Different  estimators 
lead  to  different  versions  and  implementations  of  the  PAMF  detector  [6],  [12]. 
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IV.  CG-MF  and  CG-PMF 


In  this  section,  we  discuss  alternative  implementations  of  the  MF  and  PMF  via  the  CG  algo¬ 
rithm.  The  resulting  detectors  are  referred  to  as  the  CG-MF  and  CG-PMF  detectors,  respectively, 
for  brevity.  We  start  from  the  CG-MF,  which  also  sets  the  basis  for  the  CG-PMF.  The  latter, 
by  assuming  that  the  disturbance  c(n)  is  temporally  stationary,  is  a  computationally  simplified 
version  of  the  CG-MF.  The  fink  between  the  PMF  and  CG  as  developed  in  the  sequel  reveals 
the  PMF  as  a  reduced-dimensional  subspace  detector.  In  this  section,  we  assume  knowledge 
of  the  covariance  matrix  of  the  disturbance  signal.  An  adaptive  versions  of  the  CG-PMF  (i.e., 
CG-PAMF)  is  discussed  in  Section  V. 

A.  Conjugate-Gradient  MF 

The  MF  detector,  as  shown  in  Section  III,  can  be  derived  from  a  time-varying  linear  prediction 
process.  Specifically,  consider  the  problem  of  linearly  predicting  the  n-th  sample  x(n)  under  H0 
from  all  prior  received  samples  x(n  —  1),  x(n  —  2), ... ,  x(l)  (cf.  (9)) 

x(n)  =  BH(n)y(n)  +  e(n)  (12) 

where  B(n)  =  [A^(l),  A* (2),  •  ■  ■  ,A*(n-l)]H  =  [B^n),  B2(n),  •  •  •  ,Bj(n)]  G  Cj^xj 

denotes  the  (n  —  l)-st  order  time-varying  multichannel  linear  prediction  filter,  and  y(n)  = 
[yn(l),yn(2),  ■  ■  ■  ,yn(J(n  -  1)]T  =  [xT(n  -  l),xT(n  -  2),---  ,xT(l)]T  contains  all  n  -  1 
previously  received  data  vectors.  It  is  noted  that  the  above  time-varying  linear  predictor  grows 
in  its  filter  order  or  size  with  n.  The  multichannel  linear  predictor  is  equivalent  to  J  linear 
predictors: 

Xj(n)  =  Bf  (n)y(n)  +  e^n),  j  =  1, 2,  ■  ■  ■  ,  J  (13) 

where  B j(n)  is  a  J(n  —  1) -dimensional  vector  which  contains  the  cross-channel  correlation 
information  associated  with  the  j-th  channel.  The  optimum  linear  predictor  can  be  obtained  by 
solving  the  Wiener- Hopf  equations: 

Ry(n)Bj(n)  =  Rj(n),  j  —  1,2,  ■  ■  ■  ,  J  (14) 

where  Ry(n)  =  E[y(n)yH (n)]  G  <CJ^n~1ixJ^n~v>  and  Rj(n)  =  E[y(n)x*(n)\  G  CJ(n-1)xl. 

Again,  note  that  the  size  of  the  Wiener-Hopf  equation  grows  with  n. 
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To  obtain  a  temporally  whitened  sequence  e{ri)  for  MF  detection  (cf.  (6)),  the  above  linear 
prediction  process  has  to  be  performed  multiple  times,  starting  from  n  —  2  to  n  —  N.  For  each 
n,  we  need  to  solve  a  Wiener- Flopf  equation  of  the  form  (14).  While  there  are  various  solvers  to 
the  linear  Wiener-Flopf  equation,  we  consider  using  the  conjugate  gradient  (CG)  method,  which 
has  several  properties  such  as  fast  convergence,  a  direct  link  to  the  Krylov  subspace  [7],  and 
a  built  in  model  order  selection  capability.  Additional  remarks  on  such  aspects  are  provided 
shortly. 

The  recursive  procedure  involved  for  the  determination  of  the  linear  predictors  is  described 
as  follows  (also  see  (9)). 
for  n  =  2  to  N  do 
for  j  —  1  to  ./  do 

Initialization.  Initialize  the  conjugate-direction  vector  D0y(n),  gradient  vector  7 1;-(n) 
and  initial  solution  B0j(n): 


Di  j(n)  =  =  R -jin) 

B0y  (n)  =  0. 


(15) 


(16) 


(17) 


Update  the  solution  Bfcj: 


B  kj(n)  =  B  k-ij{n)  +  akJ(n)DkJ(n). 


(18) 


Update  the  gradient  vector  7fc+1  ■: 


7fc+ij  (n)  =  7  k,j(n)  +  akij(n)Ry(n)DkJ(n). 


(19) 


Update  the  conjugate-direction  vector  Dfc+1  j-: 


(20) 


end  for 
end  for 
end  for 
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(n-i)y ,j  ( n )  !  : 

A 


Fig.  1.  Time-varying  linear  prediction  in  the  conjugate-gradient  MF  detector. 


Let  B(n)  be  the  multichannel  linear  predictor  formed  from  Bfcj  after  convergence.  Then, 

B(n)  can  be  used  to  whiten  x(n)  to  produce  a  temporally  whitened  sequence  e(n).  The  spatial 
covariance  matrix  Q(n)  of  e(n)  is  given  by  (cf.  (15)) 

Q(n)  =  E[e(n)£H  (n)] 

=  R x(n)  -  BH(n)Ryx(n)  (21) 

where  Rx(n)  =  E[x(n)x.H (n)]  E  CJxJ,  and  Ryx(n)  =  E[y(n)x.H (n)}  E  CJ(n_1)xJ,  which  is 
used  for  further  spatial  whitening  [6]. 

Fig.  1  depicts  the  CG-MF  detector  that  produces  the  n-th  sample  of  the  temporally  whitened 
sequence  £j(n)  for  the  j-th  channel,  where  D/^n)  =  [Dij(n),  D2y(n),  •  •  •  ,  Dfcj(n)]  is  the 
conjugate-direction  matrix.  CG  iterations  lead  to  a  set  of  linearly  independent  vectors  D]  j(n). . . . ,  T)k  j(n) 
that  are  conjugate  orthogonal,  i.e. 

D^(n)Ry(n)Dfi(n)  =  °>  k7 kl-  (22) 

The  output  of  the  /c-th  iteration  is  given  by 

k 

R ^  ^  dm,j(fi)Dmj(Tt)  (23) 

m=  1 

which  is  a  vector  in  the  /.--dimensional  vector  space  spanned  by  the  conjugate-direction  vectors 
{D mj(n),m  =  1,  2,  •  •  ■  The  iterative  procedure  for  the  prediction  of  the  n-th  sample  Xj(n), 
which  involves  a  J(n  —  l)-st  order  linear  predictor,  converges  after  at  most  J(n  —  1)  iterations. 

The  final  solution  Bj(n)  lies  in  a  J(n  —  1) -dimensional  vector  space. 

B.  Conjugate-Gradient  PMF  with  Known  AR  Model  Order 

If  the  disturbance  signal  can  be  approximated  as  a  temporally  wide-sense  stationary  (WSS) 
multichannel  AR  process,  the  linear  prediction  problem  of  the  previous  subsection  can  be 
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Fig.  2.  Time-invariant  linear  prediction  in  the  conjugate-gradient  PMF  detector. 


significantly  simplified.  Specifically,  suppose  the  disturbance  is  an  AR(P)  process  with  model 
order  P.  In  this  case,  the  optimum  linear  predictor  for  the  n-th  sample  x(n)  requires  only  P 
most  recently  received  samples  (as  opposed  to  all  past  samples)  and  the  prediction  filter  is 
time-invariant  with  a  fixed  size  (as  opposed  to  time- varying  with  a  growing  size)  [13]: 

x(n)  =  BHy  P(n)  +  eP{n)  (24) 

where  the  fixed  P-th  order  linear  predictor  B  =  [AH(1),AH(2),  ■  ■  ■  ,  A H(P)]H  —  [Bi,  B2,  ■  ■  ■  ,  Bj] 
G  CJPxJ  is  composed  of  the  AR  coefficient  matrices  {A  H{p)}  (cf.  (9)),  yp(n)  =  [yn(l) ,  yn(2) , 

■  ■  ■  ,yn(JP)]T  —  [xr(n  —  l),xr(n  —  2),  ■  ■  ■  ,  xr(n  —  P)]T  denotes  the  regression  data  vector, 
and  n  >  P.  Again,  it  is  convenient  to  express  the  above  multichannel  linear  predictor  as  J  scalar 
linear  predictors: 

Xj(n)  =  BfyP(n)  +  £pj(n),  i  =  1,2,  -  -  ,  J.  (25) 

The  structure  of  temporal  whitening  via  linear  prediction  for  the  PMF  detector  is  shown  in  Fig.  2. 

The  solution  to  the  scalar  linear  prediction  problem  can  be  obtained  by  solving  the  following 
Wiener-Ftopf  equation 

RyBj  =  Rj,  j  —  1,  2,  •  ■  •  ,  J  (26) 

where  Ry  =  E[yp(n)yp(n)\  G  CJPxJP  and  R7-  =  E[y p(n)x*(n)]  G  CJPxl.  It  should  be  noted 
that  unlike  the  MF  detector,  the  above  Wiener-Hopf  is  time-invariant,  has  a  fixed  size,  and  needs 
to  be  solved  only  once.  The  resulting  solution  Bj  can  be  used  to  whiten  the  entire  received  signal 
x(n)  for  n  >  P.  The  CG  algorithm  can  also  be  applied  to  solve  (26),  and  the  resulting  detector 
is  referred  to  as  the  CG-PMF  detector.  Since  only  one  fixed-sized  Wiener-Hopf  equation  needs 
to  be  solved,  the  CG-PMF  detector  is  computationally  much  simpler.  Specifically,  the  outer  loop 
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for  varying  n  as  discussed  in  Section  IV-A  vanishes,  and  only  the  conjugate-gradient  processing 
with  n  =  P  +  1  is  needed. 

Remark:  The  iterative  procedure  of  CG  converges  after  at  most  JP  iterations  for  the  CG-PMF. 
As  a  result,  the  final  solution  B,  lies  in  a  .//''-dimensional  vector  space  spanned  by  the  conjugate 
direction  vectors  Dfcj,  k  —  1,2 , ,JP,  or  equivalently,  the  .//''-dimensional  Krylov  subspace 
[7]: 

/C(Rj,  Ry,  JP)  —  span  {Rj,  RyRj,  •  •  •  ,  Ry  p_  1 R, }  .  (27) 

This  shows  that  the  PMF  is  a  reduced  JP-dimensional  solution,  as  opposed  to  the  full  JN- 
dimensional  MF  detector.  The  same  conclusion  applies  to  the  adaptive  version  CG-PAMF  detector 
discussed  in  Section  V. 


C.  Model  Order  Selection  by  CG 

In  practice,  the  AR  model  order  P  of  the  disturbance  is  often  unknown  and  has  to  be  estimated. 
A  practical  approach  is  to  choose  an  upper  bound  P  for  P,  and  use  the  CG  algorithm  to  solve 
the  following  Wiener-Flopf  equation 

Rf)Bf)  =  Rf),  j  =  1,  2,  •  •  •  ,J  (28) 


where  Ry  ’  —  E[yp(n)yp(n)]  G  cJPxJP  and  R^  =  E[yp(n)x*(n)\  G  CJPxl.  The  CG 

iterative  procedure  will  converge  after  at  most  JP  iterations  with  B^P  —  [Bp,0 /(p_p)Xl]T- 
However,  with  a  loosely  determined  upper  bound  P,  it  is  often  necessary  for  the  sake  of  reducing 
computational  complexity  to  stop  the  CG  iterations  before  it  reaches  the  maximum  number  of 
iterations.  In  this  section,  we  propose  a  model  order  selection  method  for  use  with  the  CG 
algorithm,  which  is  based  on  the  following  result. 

Lemma  1:  Suppose  the  disturbance  in  (1)  is  a  J-channel  AR(  P)  process.  Let  B  lP  G  CJPxl 
be  the  solution  to  (28)  obtained  by  CG  at  the  fc-th  iteration,  where  k  =  Jp  and  p  <  P.  Let 
By"'  G  CJpx  1  be  the  solution  to  R^By^  =  R,y"j.  Then  we  have 


—j 

- 


B[{J  =  WfcjBf ,  when  p  =  P 


(29) 


where  W/;:J  =  D|PiD'7,  D[P}  =  •  •  •  .  D is  the  conjugate-direction  matrix,  and 


D kj  is  a  k  x  k  matrix  composed  of  the  first  k  rows  of  Dfc:.y  =  D :  :I .  D2j,  •  •  •  ,  D/,j]  with 


D 


k,j 


D 


Rf'D  Z> 

(P)HTl{P)1 ^(P) 


k,j 


R-y  Dfc,j 


(30) 
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Proof:  See  Appendix  I.  ■ 

From  (29),  when  p  =  P,  the  JP  x  JP  matrix  W jpj  transforms  By  in  /C(Rr  Ry,  JP ),  which 
is  generated  by  CG-PMF  with  a  known  P,  to  B jP^  •  in  /C(Rj,  RyP\  JP),  which  is  generated  by 
CG-PMF  with  an  unknown  P.  So  the  PMF  AR  coefficient  vector  B;  is  completely  determined 
by  the  truncated  solution  B ^  •  of  CG  with  an  unknown  P  after  JP  iterations. 

We  now  explain  how  to  use  Lemma  1  for  AR  model  order  selection  in  CG-PMF.  Define  the 
residue  vector 

ekj  =  b£}  -  DgD^.Bf ,  for  k  =  Jp.  (31) 


According  to  (29),  —  0  when  k  —  JP,  so  the  norm  of  ekj  can  be  used  for  model  order 

selection.  However,  since  the  Wiener  solution  B^  is  practically  unknown,  ekj  cannot  be  directly 
computed  from  (31).  We  propose  an  approach  to  replace  B*p)  in  (31)  by  the  truncated  solution 
composed  of  the  first  k  elements  of  B[P,  which  can  be  considered  as  an  approximation  of 

[BJ,0^_P)X1F 


ek,j 


B<p)  _  t~)^Dp  B^ 

ak,j  Uk,j  Uk,jak,j 


(32) 


where  B jfj  contains  the  hrst  k  —  Jp  elements  of  B^P  .  Our  CG-based  model  order  selection 
procedure  is  summarized  as  follows. 

•  Step  1:  Select  an  upper  bound  P  for  the  model  order.  One  such  an  upper  bound  suggested 
in  [6]  for  STAP  detection  is 


P  —  max 


3  y/N 
J 


(33) 


where  |_-J  rounds  a  real-valued  number  towards  zero. 

•  Step  2:  Use  the  CG  algorithm  to  solve  the  Wiener-Hopf  equation  (28). 

-  Step  2.1:  Following  every  J  iterations  of  the  CG  algorithm,  compute  the  average  residue 
over  J  channels: 

1  . J 

=  j  "^2  llefejll  )  k  —  J,  2J, ...  (34) 

J  j= i 

-  Step  2.2:  If  e2Jp  is  smaller  than  a  specified  tolerance  level,  then  stop  the  CG  iteration, 
and  the  estimated  AR  model  order  is  P  =  p. 

The  advantage  of  the  above  CG-based  model  order  selection  method  is  that  it  does  not  require 
full  iterations  of  the  CG  algorithm  and  is  efficient.  The  complexity  of  the  CG  algorithm  with  full 
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iterations  is  in  the  same  order  as  that  of  computing  the  inverse  of  Ry'1,  which  is  0(J3P3),  while 
the  complexity  of  using  the  CG-based  order  selection  method,  is  0{J3PP2).  This  is  because  the 
latter  only  requires  JP  iterations  to  determine  the  model  order,  and  the  additional  complexity 
in  each  J  iterations  for  (32)  is  the  complexity  of  two  matrix- vector  multiplications,  which  is 
2 (JP)2-  So  the  total  complexity  is  0(JP(JP)2  +  2 P(JP)2)  «  0(J3PP2).  Next,  we  compare 
the  computational  complexity  of  the  CG-PMF  with  an  unknown  P  with  the  complexity  of  the 
eigencanceler  [14],  which  is  a  standard  eigen-subspace  detector.  The  eigencanceler  method,  has 
a  complexity  of  9J3N3  by  using  the  symmetric  QR  algorithm  to  obtain  the  eigen-subspace  and 
its  corresponding  eigen- values  [7].  Since  P  <  N,  and  generally  P  <C  N,  the  complexity  of 
CG-PMF  is  much  lower  than  eigencanceler. 

D.  Convergence  in  Airborne  Radar  Applications 

One  important  property  of  the  CG  algorithm  is  its  fast  convergence.  In  general,  it  takes  no  more 
than  JP  iterations  to  solve  the  linear  equation  (26)  [7].  Even  faster  convergence  is  possible  if  the 
covariance  matrix  of  the  disturbance  has  some  specific  structure.  In  particular,  if  the  covariance 
matrix  is  a  rank-rc  correction  of  an  identity  matrix: 

Ry  =  R,  +  ct*I  (35) 

where  R,  is  a  rank-rc  positive  semi-definite  matrix,  then  the  CG  algorithm  converges  in  at  most 
rc  +  1  iterations  [7], 

In  airborne  radar  applications,  the  disturbance  covariance  matrix  often  consists  of  two  com¬ 
ponents,  namely  a  low-rank  R,  due  to  the  clutter  and  jamming  and  a  scaled  identity  a2 I  due 
to  the  white  noise,  where  a2  denotes  the  noise  variance.  The  rank  rc  is  typically  much  smaller 
than  the  joint  spatio-temporal  dimension  JN.  Specifically,  if  the  disturbance  is  primarily  due  to 
ground  clutter  and  thermal  noise,  then  according  to  Brennan’s  rule  [2],  the  rank  of  the  clutter 
covariance  matrix  for  the  full-dimensional  MF  is  approximately 

rc,fuii  «  \J  +  {N  -  1)01  (36) 

where  f3  —  2 vgTr / d,  vg  is  the  platform  velocity,  Tr  is  the  pulse  repetition  period,  d  is  the  antenna 
element  spacing,  and  [■]  rounds  a  real-valued  number  towards  infinity. 
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Likewise,  we  can  approximate  the  rank  of  the  disturbance  covariance  matrix  for  the  PMF 
detector  as 

rcn\J  +  (P-l)P\.  (37) 

The  smaller  rank  rc  over  (36)  is  due  to  the  fact  that  the  disturbance  covariance  matrix  is  formed 
over  P  pulses,  which  is  sufficient  for  the  reduced-dimensional  PMF  detector  due  to  the  underlying 
AR( P)  model.  Meanwhile,  the  space-time  disturbance  covariance  matrix  for  the  full-dimensional 
MF  detector  is  formed  over  N  (the  entire  number  of)  pulses.  As  such,  the  PMF  can  benefit  more 
from  the  fast  convergence  property  of  the  CG  algorithm. 


E.  Preconditioned  Conjugate-Gradient  PMF 


In  cases  where  the  disturbance  covariance  does  not  have  a  low-rank  structure  as  in  (35), 
preconditioning  is  usually  helpful  in  improving  the  convergence  rate.  The  idea  is  based  on  the 
fact  that  the  convergence  rate  of  CG  is  determined  mainly  by  the  eigenvalue  structure  of  Ry. 
In  particular,  the  residue  between  the  Wiener  solution  and  fcth-step  conjugate  gradient  result  is 
bounded  by  [7] 

/  / —  -i  \  k 


llBfcd  “  BillRy  ^  2llB0,j  -  Bjlky  1 J  (38) 

where  ||o;||r.  =  y/cJ^R^cJ  denotes  the  Ry  norm  and  k  is  the  condition  number  of  Ry.  It  is 
clear  that  rapid  convergence  can  be  achieved  if  k  is  near  1.  In  the  following,  we  discuss  the  use 
of  preconditioning  with  the  CG-PMF.  For  simplicity,  the  resulting  detector  is  referred  to  as  the 
PCG-PMF  detector. 

Specifically,  consider  the  modified  Wiener-Hopf  equation  (cf.  (26)) 


RyBj  =  Rj  (39) 

where  Ry  =  M_^RyM“L  Bj  =  M^B;,  R;  =  and  M  is  a  Hermitian  positive-definite 

matrix  that  is  called  preconditioner  [7].  The  preconditioner  is  used  to  yield  a  better  conditioned 
Ry,  which  has  a  smaller  condition  number  than  Ry,  and  thus  a  faster  convergence  rate.  For 
PMF,  the  disturbance  covariance  matrix  is  a  block-Toeplitz  (BT)  matrix.  For  such  matrices, 
block-circulant  (BC)  preconditioners  are  often  recommended  [8],  [9].  Our  BC  preconditioner 
can  be  directly  computed  from  the  disturbance  covariance  matrix  Ry  which  has  the  following 
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block  Toeplitz  matrix  structure: 


where  R x(m) 
[15] 


where 


R 


y 


Rx(0) 

Rx(l  -  P) 


Rx(P-i) 

Rx(o) 


(40) 


P[x(rt,)x/f  (n  —  m)\  E  CJxJ.  In  particular,  the  BC  preconditioner  is  given  by 


M  = 


Adp  Adp_i 
Mi  M0 


Mj 

M2 


Mp_i  Mp_2  ■  ■  ■  Ado 


Adfc  = 


(P  -  fc)Rx(/c)  +  kRx{k  -  P ) 
P 


0  <  k  <  P. 


(41) 


It  is  noted  that,  as  shown  in  [7],  practically  M  2  does  not  need  to  be  explicitly  calculated  in 
the  PCG  algorithm.  The  PCG  algorithm  is  summarized  as  follows. 

Initialization.  Initialize  the  conjugate-direction  vector  D 1  ; ,  gradient  vector  7,  ■,  preconditioned 
vector  zi j  and  initial  solution  B0j: 


=  R; 

Dij  =  zij  =  M"^- 

Boy  =  0. 

for  k  —  1, 2,  ■  •  • ,  till  convergence  (k  <  J(P  —  1))  do 
Update  the  step  size  ak.f 

_  lk,j2k,j 

atJ 

Update  the  solution  Bfcj: 

B kj  B ^  1  .j  -f  j . 


Update  the  gradient  vector  -yk+ ,  ;  : 

Tfc+iy  =  TTy  T  OfcjRyD^j. 


(42) 

(43) 

(44) 


(45) 

(46) 
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(47) 


Update  the  preconditioned  vector  Zk+i.j' 

Zfc+ij  =  M_17fc+1J.  (48) 

Update  the  conjugate-direction  vector  Dfc+]  j 


~f~ 


'Tfe+l,jzfc+l,j 


Dfc  j . 


(49) 


end  for 

The  complexity  associated  with  the  AR  parameter  estimation  in  PCG-PMF  is  summarized 
in  TABLE  I,  where  r  is  the  number  of  iterations  needed  by  the  PCG  algorithm  to  reach 
convergence,  and  the  flop  counts  are  for  all  J  channels.  It  is  interesting  to  note  that  the  PCG-PMF 
is  computationally  very  efficient,  involving  approximately  0(rJ3P  log2  P).  The  computational 
efficiency  is  primarily  due  to  the  fast  convergence  rate  offered  by  preconditioning  and  the  use 
of  a  BC  preconditioner,  as  explained  next.  In  the  following,  we  discuss  the  complexity  of  only 
M-1,  (45)  and  (48),  since  the  other  calculations  are  obvious. 

First,  we  consider  M-1.  Since  M  is  a  block-circulant  (BC)  matrix,  the  inverse  of  M  can  be 
computed  by  using  the  Fast  Fourier  transform  (FFT)  [16] 


where 


C0 


Cp-i 


Cp_!  •••  C, 
Co  •  •  •  C2 

C  P_2  Co 


Cm 


1 

P2 


P-1 

J]  IUpfcmM^ 

k= 0 


m  =  0, 1,  •  •  •  ,P  —  1 


(50) 


(51) 


and  WPkm  —  exp(j2kmn / P) .  It  follows  that  the  computation  of  M  1  is  composed  of  P  matrix 
inversions  of  J  x  J  matrices  and  J2  FFTs  of  length  P.  Therefore,  the  total  complexity  is 
0{J2P  log2  P  +  J3P). 

Second,  we  consider  (45).  The  main  complexity  of  (45)  is  matrix-vector  multiplication  RyD/  r 
Since  Ry  is  a  JP-dimensional  BT  matrix,  the  above  matrix-vector  multiplication  consists  of  J2 
Toeplitz  matrix-vector  multiplications,  where  each  Toeplitz  matrix  is  a  P  x  P  matrix.  Each 
Toeplitz  matrix-vector  multiplication  can  be  implemented  by  the  FFT  using  0(P\og2P)  flops 
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TABLE  I 


Complexity  of  PCG-PMF 


Equation 

Flops 

Remark 

(41) 

0(J2P ) 

calculated  once 

M-1 

O {J2 P  \og2  P  +  J3P) 

calculated  once 

(45) 

0{J3P  log2  P ) 

at  fc-th  iteration 

(46) 

0{J2P) 

at  fc-th  iteration 

(47) 

0(J2P) 

at  fc-th  iteration 

(48) 

0(J3P  log2  P) 

at  fc-th  iteration 

(49) 

0(J2P) 

at  fc-th  iteration 

Total 

«  0{rJ3P  log2  P) 

for  r  iterations 

[17].  Hence,  the  complexity  of  (45)  for  each  channel  per  iteration  is  0(J2P  \og2  P).  With  J 
channels  and  r  iterations,  the  total  complexity  of  (45)  is  0(rJ3P  log2  P). 

Finally,  we  consider  (48).  Since  the  preconditioner  M  is  a  BC  matrix,  (48)  can  again  be  com¬ 
puted  by  J2  FFTs  of  length  P.  The  complexity  for  each  channel  per  iteration  is  0(J2P  log2  P), 
so  the  total  complexity  of  (48)  for  J  channels  is  0{rJ3P  log2  P ). 

Here,  we  make  a  comparison  between  the  PCG-PMF  and  CG-PMF.  Since  the  condition  number 
of  the  preconditioned  disturbance  covariance  matrix  Ry  is  generally  smaller  than  that  of  Ry, 
PCG-PMF  provides  a  faster  convergence  than  CG-PMF.  The  latter  has  a  complexity  of  0(J3P3). 

V.  Conjugate-Gradient  PAMF 

The  CG-PMF  algorithm  is  now  extended  to  the  adaptive  case  when  both  the  covariance  matrix 
and  the  AR  model  order  P  are  unknown.  The  resulting  detector  is  referred  to  as  the  CG-PAMF 
detector.  The  extension  of  CG-PMF  involves  i)  replacing  the  true  covariance  matrices  with 
estimates  obtained  from  the  target-free  training  data;  and  ii)  integrating  the  CG-based  model 
order  selection  proposed  in  Section  IV-C  with  conjugate-gradient  iterations.  The  CG-PAMF 
detector  with  order  selection  (OSCG-PAMF)  is  summarized  next. 

•  Step  1:  Estimate  the  disturbance  covariance  matrices  from  the  training  data  via  temporal 
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and  range  averaging: 


Rx(0)  •  •  •  R  X(P  -  1) 


(52) 


Rx(l  -P)  ■  •  ■  Rx(0) 


(53) 


where  the  sub-matrices  are  given  by 


i  \  —  in 

RxxM  =  jfx  S  J2  *k(l  + 


K  N-m 


(54) 


k= 1  1=1 


with  K  denoting  the  number  of  training  data  vectors  and  P  determined  by  (33). 
•  Step  2:  Use  the  CG  algorithm  to  solve 


(55) 


-  Step  2.1:  Examine  the  residual  ejp  (34)  at  each  Jp- th  (p  —  1,  2,  ■  ■  ■  ,  P)  iteration  of 
CG.  If  ejp  <  Q'oejfp-i),  where  0  <  a0  <  !  is  a  stopping  threshold,  then  exit  the  CG 
iteration,  and  set  the  AR  model  order  as  P  =  p. 


Unlike  the  original  PAMF  with  an  unknown  AR  model  order  [6],  which  has  to  run  recursively 
from  p  —  1  to  a  P  (P  <  P)  to  jointly  estimate  the  AR  coefficients  and  model  order,  OSCG- 


PAMF  does  not  contain  any  recursion.  It  only  has  to  perform  CG  with  the  disturbance  covariance 


(P) 

matrix  Ry  for  JP  iterations  to  obtain  a  model  order  estimate. 


Remark:  Several  estimators  can  be  employed  to  find  the  linear  prediction  filters  for  the  PAMF. 
The  estimator  as  represented  by  (55)  along  with  the  covariance  matrix  estimates  (52)-(54)  is 
often  called  the  multichannel  Yule- Walker  method.  Other  estimators,  such  as  the  least-squares 
estimators  [6],  solve  slightly  modified  versions  of  the  linear  equation  (55).  It  is  noted  that  in 
most  cases,  the  CG  algorithm  can  be  used  to  efficiently  solve  such  a  modified  linear  equation. 
Due  to  space  limit,  we  will  not  explore  these  alternative  CG-PAMF  detectors. 

A  similar  comparison  can  be  made  between  the  complexity  of  the  OSCG-PAMF  detector  and 
that  of  the  eigencanceler  when  the  covariance  matrix  is  unknown.  In  addition  to  the  numbers  of 
flops  as  summarized  in  Section  IV-C,  both  have  to  pay  the  extra  complexity  needed  to  estimate 
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the  covariance  matrix.  In  this  case,  the  OSCG-PAMF  requires  an  additional  complexity  of 
0(J2PNK )  as  incurred  in  (52)-(54),  whereas  the  extra  complexity  for  the  eigencanceler  is 
0(J2N2K)  that  is  used  to  estimate  a  full  (JN  x  JN)  space-time  covariance  matrix  from  K 
training  signals. 


VI.  Numerical  Results 

In  this  section,  simulation  results  are  provided  to  illustrate  the  performance  of  the  proposed 
techniques.  We  consider  simulated  data  generated  using  AR  models  and  the  KASSPER  data 
[18],  which  were  obtained  from  more  realistic  clutter  models.  The  simulation  results  presented 
below  use  20000  independent  Monte  Carlo  data  realizations  and  a  probability  of  false  alarm  of 
.Pfa  =  10-2.  The  chosen  Pf.d  may  be  considered  too  high  for  many  practical  detection  applications. 
It  is  noted  that  the  choice  is  only  to  facilitate  computer  simulation  and  reduce  simulation  time. 
The  main  observations  from  the  simulation,  including  the  convergence  of  the  CG  algorithm  in 
PAMF  detection  and  the  accuracy  of  the  estimated  P  provided  by  the  proposed  model  order 
selection  method,  are  independent  of  the  choice  of  Pfa. 

A  major  issue  that  we  like  to  illustrate  in  the  following  numerical  examples  is  the  convergence 
of  the  CG  algorithm,  with  partial  or  full  iterations,  when  the  data  covariance  matrix  is  known 
or  estimated,  and/or  when  the  AR  model  order  is  known  or  estimated.  To  this  end,  we  compare 
the  various  detectors  used  with  the  CG  algorithm,  including  the  CG-PMF  (Section  IV-B),  CG- 
PAMF  with  a  known  AR  model  order  (Section  V)  and  OSCG-PAMF  with  an  estimated  model 
order  (Section  V),  with  the  same  detectors  used  with  direct  matrix  inverse  (DMI).  For  example, 
the  DMI-PAMF  detector  involves  a  direct  inverse  of  the  estimated  covariance  matrix  in  (52) 
and  uses  it  compute  the  linear  prediction  hlter  (26).  This  DMI  approach  turns  out  to  coincide 
with  the  Yule- Walker  method  [13]  for  AR  spectral  estimation.  It  is  noted  in  [6]  that  there  are 
alternative  spectral  estimation  methods  which  may  yield  better  detection  performance  in  some 
scenarios.  These  alternatives  are  not  considered  here  since  the  focus  is  the  convergence  of  the 
CG  algorithm  in  PAMF.  In  the  following,  we  will  primarily  use,  as  a  comparison  metric,  the 
probability  of  detection  versus  the  SINR  for  a  given  probability  of  false  alarm.  The  output  SINR 
of  the  PAMF  detector  was  derived  and  extensively  studied  in  [19]. 

First,  we  examine  the  performance  of  the  two  implementations  of  the  PMF  detector  by  using 
simulated  data  with  AR  disturbances.  The  disturbance  is  an  AR(2)  process  with  J  —  4  elements 
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J=4,  N=64,  P  =1 .000000e-002,  AR(2) 
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0  5  10 

SINR  (dB) 


Fig.  3.  Probability  of  detection  versus  SINR  of  PMF  for  simulated  data(J  =  4;  N  =  64;  P  =  2) 


and  N  =  64  pulses.  Both  PMF  detectors  have  knowledge  of  the  exact  disturbance  covariance 
matrix;  however,  they  use  different  approaches  to  compute  the  linear  predictor.  Specihcally,  we 
consider  the  DMI-PMF,  which  uses  direct  matrix  inverse  to  solve  the  Wiener-FIopf  equation, 
and  the  CG-PMF  as  discussed  in  Section  IV-B  with  knowledge  of  the  AR  model  order  P.  The 
numerical  results  are  shown  in  Fig.  3.  It  is  seen  that  both  implementations  yield  an  identical 
detection  performance. 

We  next  examine  the  performance  of  the  CG-based  AR  model  order  selection  method  used  in 
the  CG-PMF  and  CG-PAMF  detectors  with  an  unknown  AR  model  order.  Two  AR  disturbance 
signals  with  J  —  4  and  N  =  63  are  considered,  and  their  model  orders  are  P  —  1  and  3, 
respectively.  We  choose  the  same  upper  bound  P  —  6  for  both  cases.  The  residual  ck  (34)  is 
computed  and  used  for  model  order  selection;  as  a  benchmark,  we  also  include  efe,  which  is 
similarly  computed  as  in  (34)  but  with  ekj  replaced  by  ekj.  Recall  that  q.  is  an  approximation 
of  Ik,  which  cannot  be  computed  in  practice  due  to  the  fact  that  the  true  Wiener-Hopf  solution  is 
unknown.  The  numerical  results  for  CG-PMF  are  shown  in  Fig.  4  and  Fig.  5,  which  correspond 
to  P  —  1  and  P  =  3,  respectively.  It  is  seen  that,  ek  has  a  sharp  decrease  at  the  JP-th  ( JP  —  4 
for  P  —  1  and  JP  —  12  for  P  —  3)  iteration  of  CG,  which  confirms  the  effectiveness  of  the 
CG-based  model  order  selection  method.  The  counterpart  model  order  selection  results  for  CG- 
PAMF  are  shown  in  Fig.  6  and  Fig.  7,  for  P  =  1  and  P  —  3,  respectively.  Unlike  the  CG-PMF 
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Fig.  4.  Residuals  ejv  and  ejp  for  model  order  selection  in  CG-PMF  (J  =  4;  N  =  64;  P  =  1;  P  =  6) 
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Fig.  5.  Residuals  ejp  and  ijp  for  model  order  selection  in  CG-PMF  (J  s*  4;  N  =  64;  P  =  3;  P  =  6) 


which  uses  the  real  disturbance  covariance  matrix,  the  sample  covariance  matrix  estimated  from 
the  training  data  (cf.  (40))  is  employed  for  model  order  selection  in  CG-PAMF.  Here  the  training 
data  size  is  set  to  K  —  32.  It  is  also  shown  that  ck  has  a  sharp  decrease  at  the  JP-th  ( JP  =  4 
for  P  —  1  and  -IP  —  VI  lor  P  =  3)  iteration  of  CG,  although  the  decrease  in  residue  is  smaller 
than  that  of  CG-PMF  due  to  estimation  error  of  the  sample  covariance  matrix. 

We  now  consider  the  convergence  of  PCG-PMF.  The  simulated  disturbance  is  an  AR(8) 
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K=32,  J=4,  N=64,  P  =6,  AR(1 ) 


Fig.  6.  Residuals  ejp  and  ejp  for  model  order  selection  in  CG-PAMF  ( K  =  32;  J  =  4;  N  =  64;  P  =  1;  P  =  6) 


K=32,  J=4,  N=64,  P  =6,  AR(3) 


Number  of  iterations  (k) 


Fig.  7.  Residuals  ejp  and  ijp  for  model  order  selection  in  CG-PAMF  (K  =  32;  J  =  4;  N  =  64;  P  =  3;  P  =  6) 


multichannel  process  with  J  —  4.  The  convergence  of  CG-PMF  and  PCG-PMF  is  shown  in 
Fig.  8.  The  condition  number  of  the  preconditioned  covariance  matrix  is  4.2,  which  is  much  less 
than  the  condition  number  of  the  original  covariance  matrix  77.1.  It  is  seen  from  Fig.  8  that  only 
5  iterations  are  needed  in  PCG-PMF  to  reach  a  relative  approximation  error  under  1%,  while  20 
iterations  are  needed  for  CG-PMF. 

Our  next  example  considers  the  adaptive  PAMF  detector,  for  which  the  disturbance  covariance 
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Fig.  8.  Convergence  of  CG-PMF  and  PCG-PMF  ( J  =  4;  P  =  8) 


matrix  is  unknown  and  the  sample  covariance  matrix  is  estimated  by  (52).  Similar  to  the  PMF 
detector,  we  compare  two  implementations  of  the  PAMF  detector,  including  the  DMI-PAMF 
and  CG-PAMF,  Flere  the  DMI-PAMF  directly  inverses  the  sample  covariance  matrix  to  get  the 
maximum-likelihood  estimation  of  AR  coefficients  [13].  The  disturbance  is  an  AR(2)  signal, 
whose  disturbance  covariance  matrix  is  estimated  from  K  =  16  target- free  training  data  vectors, 
and  the  AR  coefficients  are  estimated  based  on  the  estimated  disturbance  covariance  matrix.  The 
numerical  results  are  shown  in  Fig.  9.  It  is  observed  that  both  implementations  yield  an  identical 
detection  performance. 

The  performance  of  the  CG-PAMF  with  an  unknown  disturbance  AR  model  order  and  dis¬ 
turbance  covariance  matrix  is  considered  next.  Both  AR  model  based  data  and  KASSPER  2002 
data  set  are  employed  in  this  example.  The  KASSPER  data  set  was  generated  by  considering 
practical  airborne  radar  parameters  and  issues  found  in  a  real-world  clutter  environment  [18]. 
Specifically,  the  simulated  airborne  radar  platform  travels  at  a  speed  of  100  m/s  with  a  3°  crab 
angle.  The  radar  carrier  frequency  is  1240  MHz.  The  horizontal  11  antenna  elements  form  a 
ULA  with  a  spacing  of  0.1092m  between  adjacent  elements,  and  the  transmit  array  is  uniformly 
weighted  and  phased  to  steer  the  main  beam  to  195°.  The  pulse  repetition  frequency  is  1984 
Hz  and  a  coherent  processing  interval  contains  32  pulses.  Only  the  first  8  elements  are  used  in 
our  simulation.  We  use  the  covariance  matrix  associated  with  range  bin  200  in  the  KASSPER 


22 

Approved  for  public  release;  distribution  unlimited. 


K=16,  J=4,  N=64,  P  =1 .000000e-002,  AR(2) 


Fig.  9.  Probability  of  detection  versus  SINR  of  PAMF  for  AR  data  ( K  =  16;  J  =  4;  N  =  64;  P  =  2) 


data  set  to  generate  the  test  data  and  the  covariance  matrices  from  the  neighboring  ranges  bins 
to  generate  the  training  signals.  A  target  is  injected  into  the  test  cell  with  a  normalized  spatial 
frequency  0.1  and  a  normalized  Doppler  frequency  0.35. 

The  numerical  results  are  shown  in  Fig.  10  for  the  AR  model  based  data  and,  respectively, 
Fig.  11  for  the  KASSPER  2000  data,  where  OSCG-PAMF  (unknown  P )  represents  the  CG- 
PAMF  detector  with  the  CG-based  model  order  selection  method  which  employs  a  model  order 
upper  bound  P  calculated  by  (33).  In  Fig.  11,  we  also  include  for  comparison  the  joint  domain 
localized  (JDL)  detector  [20],  a  popular  reduced-dimensional  STAP  solution  in  scenarios  of 
limited  training.  The  IDL  is  implemented  by  using  3  beams  and  3  Doppler  bins  for  adaptivity. 
It  is  seen  that  the  performance  of  the  OSCG-PAMF  is  nearly  identical  to  that  of  CG-PAMF  with 
known  P  (AR  data)  or  a  pre-selected  P  —  2  (KASSPER  data).  For  the  case  of  AR  data,  we 
noticed  that  only  one  model  order  selection  error  (P  /  P)  occurred  out  of  20000  simulations. 
Moreover,  using  the  relevant  parameters  of  the  KASSPER  data,  we  have  /3  —  2 vgTr/d  —  0.923. 
It  follows  that  for  J  —  8  elements,  the  maximum  number  of  conjugate-gradient  iterations  needed 
by  the  CG-PAMF  for  a  given  model  order  p  is  estimated  to  be  rcp  +  1  =  [8.077  +  0.923/;] .  For 
example,  the  maximum  numbers  of  CG  iterations  for  p  =  2  is  about  10  due  to  the  low -rank 
structure  of  the  clutter,  whereas  without  such  a  structure,  it  would  require  pj  —  16  iterations  for 
the  CG  to  converge.  It  is  also  seen  in  Fig.  1 1  that  the  PAMF  detectors  outperform  the  JDF-AMF 
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K=32,  J=4,  N=64,  P  =1 .000000e-002,  P  =6 

f  max 


Fig.  10.  Probability  of  detection  versus  SINR  for  AR  data  ( K  =  32;  J  =  4;  N  =  64;  P  =  2;  P  =  6) 


K=32,  J=8,  N=32,  Pf=1 .000000e-002,  Pmax=2 


Fig.  11.  Probability  of  detection  versus  SINR  for  KASSPER  2002  data  (K  =  32;  J  =  8;  N  =  32;  P  =  2;  P  =  2) 


detector.  The  JDL-AML  experiences  a  loss  of  about  4  dB  compared  with  the  MF. 

Finally,  we  compare  the  complexity  in  terms  of  the  number  of  flops  required  by  the  CG  and 
DMI  implementations.  The  flops  required  by  the  CG-PAMF  and  DMI-PAMF  versus  the  AR 
model  order  p  are  shown  in  Fig.  12.  For  the  DMI-PAMF,  the  QR  decomposition  is  adopted  to 
get  the  J-channel  AR  coefficients.  It  is  seen  that  the  complexity  of  the  CG-PAMF  is  lower  than 
that  of  the  DMI-PAMF. 
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Fig.  12.  Computational  complexity  of  CG-PAMF  and  DM1-PAMF  versus  AR  model  order  p  (K  =  32;  J  =  8;  N  =  32). 


VII.  Conclusions 

The  conjugate-gradient  (CG)  algorithm  was  employed  to  solve  the  linear  prediction  problem 
underlying  the  parametric  matched  filter  (PMF)  and  parametric  adaptive  matched  filter  (PAMF) 
detectors.  It  is  shown  that  the  CG  algorithm  leads  to  not  only  new  efficient  implementations, 
but  also  new  insights  of  these  parametric  detectors  as  reduced-dimensional  subspace  detectors. 
In  particular,  the  linear  prediction  filter  and  whitening  filter  of  the  PMF  and  PAMF  detectors  are 
within  the  Krylov  subspace  of  dimension  JP,  and  these  detectors  are  reduced  .//^-dimensional 
subspace  detectors,  where  J  and  P  are  the  number  of  channels  and  AR  model  order,  respec¬ 
tively.  We  examined  the  convergence  rate  of  the  CG  parametric  detectors.  In  airborne  radar 
applications,  the  special  low-rank  structure  of  the  disturbance  covariance  matrix  implies  that  a 
rapid  convergence  is  possible,  whereby  convergence  can  be  achieved  without  completing  a  full 
round  of  CG  iterations.  Even  for  disturbance  covariance  matrices  that  do  not  have  the  low -rank 
structure,  preconditioning  methods  can  be  used  to  speed  up  the  convergence  rate.  In  general,  the 
CG  parametric  detectors  are  more  efficient  than  their  counterparts  implemented  in  conventional 
approaches.  We  also  presented  a  new  CG-based  AR  model  order  selection  method,  which  is 
naturally  integrated  with  the  CG  iterations.  The  proposed  techniques  are  illustrated  by  using 
both  KASSPER  and  other  simulated  data. 

Finally,  we  note  that  the  CG  algorithm  bears  some  similarity  to  a  vector  space  approach  [21] 
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to  solving  the  multi-dimensional  Yule- Walker  equation  for  an  arbitrary  region  of  support.  Both 
involve  the  use  of  conjugate  orthogonal  basis  vectors.  A  future  subject  would  be  to  investigate 
the  relation  of  the  two  approaches  and  explore  the  application  of  the  CG  algorithm  for  multi¬ 
dimensional  and  multichannel  applications. 


Appendix  I 
Proof  of  Lemma  1 

It  is  well  known  that  the  conjugate  direction  vectors  obtained  by  the  CG  algorithm  solving 
the  Wiener- Hopf  equation  (28)  span  the  Krylov  subspace  [7]: 


/C(Rf Rf),  k)  =  span  {d<  j,  D^,  ■  •  •  ,  Djg} 


(56) 


Furthermore,  the  truncated  solution  obtained  at  the  fc-th  iteration  minimizes  the  Ry '  -norm 


of  the  approximation  error  over  all  vectors  on  /C(Rj,  RyP\  k)  [22],  i.e. 
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Therefore,  the  truncated  solution  obtained  at  the  fc-th  iteration  Bj,.  .  is  the  Ry  -orthogonal  projec¬ 
tion  of  the  Wiener  solution  B*P^  to  the  subspace  /C(Rj,  RyP),  k),  and  cxkj  — 
contains  the  coordinate  values  of  conjugate-direction  vectors  |d^, 
given  by 


d-G 

C  ’ 


,  D^J,  which  are 


a(P)  = 
kj 


(58) 


With  the  definition  of  by  (30),  we  can  write  after  JP  iterations 

Cjpj  =  OjPjBf' 


(59) 


a-JPJ  - 


^ JP,j  ^JP,  d 


where  t)JPJ  =  [DUr  D 2J,  •  ■  ■  ,  DjPJ].  Recalling  Bp  =  [Bj,  0j(J5_P)xl] 

B, 


,  we  have 


0 


(60) 


where  DjPj  G  CJPxJP  is  the  upper  JP  x  JP  block  matrix  of  Djpy,  and  DjpjCi  contains  the 
lower  block  of  DjPj.  Then  Bj  is  given  by 


B, 


^jpjaJP,j- 


(61) 
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The  intermediate  solution  obtained  at  the  JP-th  CG  iteration  is 

jp 

bj%  =  E  <wD! 3  =  <62> 

m=l 

It  follows  from  (61)  and  (62)  that  B  {fp  ■  and  B j  are  related  by 

=  Dg  WJPjBj  (63) 

where  Wj Pj  —  jD,JVj  E  CJT‘/  JP,  which  completes  the  proof. 
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